{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/W2D1-postcourse-bugfix/tutorials/W2D2_LinearSystems/student/W2D2_Tutorial2.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "# Neuromatch Academy 2020, Week 2, Day 2, Tutorial 2\n",
    "\n",
    "# Markov Processes\n",
    "\n",
    "**Content Creators**: Bing Wen Brunton, Ellie Stradquist\n",
    "\n",
    "**Content Reviewers**: Norma Kuhn, Karolina Stosio, John Butler, Matthew Krause, Ella Batty, Richard Gao, Michael Waskom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Tutorial Objectives\n",
    "\n",
    "In this tutorial, we will look at the dynamical systems introduced in the first tutorial through a different lens. \n",
    "\n",
    "In Tutorial 1, we studied dynamical systems as a deterministic process. For Tutorial 2, we will look at **probabilistic** dynamical systems. You may sometimes hear these systems called _stochastic_. In a probabilistic process, elements of randomness are involved. Every time you observe some probabilistic dynamical system, started from the same initial conditions, the outcome will likely be different. Put another way, dynamical systems that involve probability will incorporate random variations in their behavior. \n",
    "\n",
    "For some probabilistic dynamical systems, the differential equations express a relationship between $\\dot{x}$ and $x$ at every time $t$, so that the direction of $x$ at _every_ time depends entirely on the value of $x$. Said a different way, knowledge of the value of the state variables $x$ at time t is _all_ the information needed to determine $\\dot{x}$ and therefore $x$ at the next time.\n",
    "\n",
    "This property --- that the present state entirely determines the transition to the next state --- is what defines a **Markov process** and systems obeying this property can be described as **Markovian**.\n",
    "\n",
    "The goal of Tutorial 2 is to consider this type of Markov process in a simple example where the state transitions are probabilistic. In particular, we will:\n",
    "\n",
    "* Understand Markov processes and history dependence.\n",
    "* Explore the behavior of a two-state telegraph process and understand how its equilibrium distribution is dependent on its parameters.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "#@title Figure settings\n",
    "import ipywidgets as widgets       # interactive display\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "#@title Helper Functions\n",
    "\n",
    "def plot_switch_simulation(t, x):\n",
    "  fig = plt.figure()\n",
    "  plt.plot(t, x)\n",
    "  plt.title('State-switch simulation')\n",
    "  plt.xlabel('Time')\n",
    "  plt.xlim((0, 300)) # zoom in time\n",
    "  plt.ylabel('State of ion channel 0/1', labelpad=-60)\n",
    "  plt.yticks([0, 1], ['Closed (0)', 'Open (1)'])\n",
    "  plt.show()\n",
    "  return\n",
    "\n",
    "def plot_interswitch_interval_histogram(inter_switch_intervals):\n",
    "  fig = plt.figure()\n",
    "  plt.hist(inter_switch_intervals)\n",
    "  plt.title('Inter-switch Intervals Distribution')\n",
    "  plt.ylabel('Interval Count')\n",
    "  plt.xlabel('time')\n",
    "  plt.show()\n",
    "\n",
    "def plot_state_probabilities(time, states):\n",
    "  fig = plt.figure()\n",
    "  plt.plot(time, states[:,0], label='Closed to open')\n",
    "  plt.plot(time, states[:,1], label='Open to closed')\n",
    "  plt.legend()\n",
    "  plt.xlabel('time')\n",
    "  plt.ylabel('prob(open OR closed)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 1:  Telegraph Process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 519
    },
    "colab_type": "code",
    "outputId": "d3b9f2a8-8e88-4f34-8bab-9deca3650ec2"
   },
   "outputs": [],
   "source": [
    "#@title Video 1: Markov Process\n",
    "# Insert the ID of the corresponding youtube video\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"xZO6GbU48ns\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Let's consider a Markov process with two states, where switches between each two states are probabilistic (known as a telegraph process). To be concrete, let's say we are modeling an **ion channel in a neuron that can be in one of two states: Closed (0) or Open (1)**. \n",
    "\n",
    "If the ion channel is Closed, it may transition to the Open state with probability $P(0 \\rightarrow 1 | x = 0) = \\mu_{c2o}$. Likewise, If the ion channel is Open, it transitions to Closed with probability $P(1 \\rightarrow 0 | x=1) = \\mu_{o2c}$.\n",
    "\n",
    "We simulate the process of changing states as a **Poisson process**. The Poisson process is a way to model discrete events where the average time between event occurrences is known but the exact time of some event is not known. Importantly, the Poisson process dictates the following points: \n",
    "1. The probability of some event occurring is _independent from all other events_.\n",
    "2. The average rate of events within a given time period is constant.\n",
    "3. Two events cannot occur at the same moment. Our ion channel can either be in an open or closed state, but not both simultaneously. \n",
    "\n",
    "In the simulation below, we will use the Poisson process to model the state of our ion channel at all points $t$ within the total simulation time $T$. \n",
    "\n",
    "As we simulate the state change process, we also track at which times thoughout the simulation the state makes a switch. We can use those times to measure the distribution of the time _intervals_ between state switches."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "**Run the cell below** to show the state-change simulation process. Note that a random seed was set in the code block, so re-running the code will produce the same plot. Commenting out that line will produce a different simulation each run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "outputId": "5fad8266-7913-4c70-a085-0087533a83ac"
   },
   "outputs": [],
   "source": [
    "# @title State-change simulation process\n",
    "\n",
    "# parameters\n",
    "T = 5000     # total Time duration\n",
    "dt = 0.001   # timestep of our simulation\n",
    "\n",
    "# simulate state of our ion channel in time\n",
    "# the two parameters that govern transitions are\n",
    "# c2o: closed to open rate\n",
    "# o2c: open to closed rate\n",
    "def ion_channel_opening(c2o, o2c, T, dt):\n",
    "    # initialize variables\n",
    "    t = np.arange(0, T, dt)\n",
    "    x = np.zeros_like(t)\n",
    "    switch_times = []\n",
    "\n",
    "    # assume we always start in Closed state\n",
    "    x[0] = 0\n",
    "\n",
    "    # generate a bunch of random uniformly distributed numbers\n",
    "    # between zero and unity: [0, 1),\n",
    "    # one for each dt in our simulation.\n",
    "    # we will use these random numbers to model the\n",
    "    # closed/open transitions\n",
    "    myrand = np.random.random_sample(size=len(t))\n",
    "\n",
    "\n",
    "    # walk through time steps of the simulation\n",
    "    for k in range(len(t)-1):\n",
    "        # switching between closed/open states are\n",
    "        # Poisson processes\n",
    "        if x[k] == 0 and myrand[k] < c2o*dt: # remember to scale by dt!\n",
    "            x[k+1:] = 1\n",
    "            switch_times.append(k*dt)\n",
    "        elif x[k] == 1 and myrand[k] < o2c*dt:\n",
    "            x[k+1:] = 0\n",
    "            switch_times.append(k*dt)\n",
    "\n",
    "    return t, x, switch_times\n",
    "\n",
    "\n",
    "c2o = 0.02\n",
    "o2c = 0.1\n",
    "np.random.seed(0) # set random seed\n",
    "t, x, switch_times = ion_channel_opening(c2o, o2c, T, .1)\n",
    "plot_switch_simulation(t,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 1 (2A): Computing intervals between switches\n",
    "We now have `switch_times`, which is a list consisting of times when the state switched. Using this, calculate the time intervals between each state switch and store these in a list called `inter_switch_intervals`.\n",
    "\n",
    "We will then plot the distribution of these intervals. How would you describe the shape of the distribution?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "##############################################################################\n",
    "## TODO: Insert your code here to calculate between-state-switch intervals,\n",
    "##       and uncomment the last line to plot the histogram\n",
    "##############################################################################\n",
    "\n",
    "# hint: see np.diff()\n",
    "# inter_switch_intervals = ...\n",
    "\n",
    "\n",
    "# plot_interswitch_interval_histogram(inter_switch_intervals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 433
    },
    "colab_type": "text",
    "outputId": "407f2b47-bc20-4723-9f97-f5ae3d2d98fb"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D2_LinearSystems/solutions/W2D2_Tutorial2_Solution_30701e58.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=416 src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W2D2_LinearSystems/static/W2D2_Tutorial2_Solution_30701e58_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "We can also generate a bar graph to visualize the distribution of the number of time-steps spent in each of the two possible system states during the simulation. **Run the cell below** to visualize the distribution. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "outputId": "1b4f1d75-b59b-4524-e167-ba6de9a14d7f"
   },
   "outputs": [],
   "source": [
    "# @title Distribution of time spent in each state.\n",
    "\n",
    "states = ['Closed', 'Open']\n",
    "(unique, counts) = np.unique(x, return_counts=True)\n",
    "plt.bar(states, counts)\n",
    "plt.ylabel('Number of time steps')\n",
    "plt.xlabel('State of ion channel');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "<!-- Though the system started initially in the Closed ($x=0$) state, over time, it settles into a equilibrium distribution where we can predict on what fraction of time it is Open as a function of the $\\mu$ parameters. \n",
    "\n",
    "Before we continue exploring these distributions further, let's first take a look at the this fraction of Open states as a cumulative mean of the state $x$: -->\n",
    "\n",
    "Even though the state is _discrete_--the ion channel can only be either Closed or Open--we can still look at the **mean state** of the system, averaged over some window of time. \n",
    "\n",
    "Since we've coded Closed as $x=0$ and Open as $x=1$, conveniently, the mean of $x$ over some window of time has the interpretation of **fraction of time channel is Open**.\n",
    "\n",
    "Let's also take a look at the fraction of Open states as a cumulative mean of the state $x$. The cumulative mean tells us the average number of state-changes that the system will have undergone after a certain amount of time. **Run the cell below**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 431
    },
    "colab_type": "code",
    "outputId": "cab9f14b-ae8d-428d-cc57-aa737f99e962"
   },
   "outputs": [],
   "source": [
    "# @title Cumulative mean of state\n",
    "plt.plot(t, np.cumsum(x) / np.arange(1, len(t)+1))\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('Cumulative mean of state');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Notice in the plot above that, although the channel started in the Closed ($x=0$) state, gradually adopted some mean value after some time. This mean value is related to the transition probabilities $\\mu_{c2o}$\n",
    "and $\\mu_{o2c}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: Varying transition probability values & T\n",
    "\n",
    "Using the interactive demo below, explore the state-switch simulation for different transition probability values of states $\\mu_{c2o}$ and $\\mu_{o2c}$. Also, try different values for total simulation time length *T*. \n",
    "\n",
    "Does the general shape of the inter-switch interval distribution change or does it stay relatively the same? How does the bar graph of system states change based on these values?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 526,
     "referenced_widgets": [
      "0cb53bd8ab4c4371b28166f3b89f0c68",
      "fc03e6bf6c4841f8ac0f20746d302e68",
      "5aa7136ed3af4359b594468ca0b9a99b",
      "37cbcc95605f44659f99467a9e9fcbfd",
      "bec5b78e22a342e0ba733362d6d5561b",
      "e32b1207067c40dabead946520475a43",
      "5a9ee70b015d46ac9219bf7b50e17307",
      "11ee88f4d3df4d3fb61b3d33254b89f4",
      "afce7d7951ed4e629fac35b4c6b40d89",
      "98ecccfbe7eb469babdc7437f5a8b7a2",
      "485e8f8f88014f838f5305984ff67c41",
      "f4f5aa61772d4cf583afb59c66ca35e6",
      "2810d78a37dd46c686c7ee1397d28697"
     ]
    },
    "colab_type": "code",
    "outputId": "8d82d162-6d77-432c-a3a8-74dd74ffa10c"
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "\n",
    "#@markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "@widgets.interact\n",
    "def plot_inter_switch_intervals(c2o = (0,1, .01), o2c = (0, 1, .01), T=(1000,10000, 1000)):\n",
    "\n",
    "  t, x, switch_times = ion_channel_opening(c2o, o2c, T, .1)\n",
    "\n",
    "  inter_switch_intervals = np.diff(switch_times)\n",
    "\n",
    "  #plot inter-switch intervals\n",
    "  plt.hist(inter_switch_intervals)\n",
    "  plt.title('Inter-switch Intervals Distribution')\n",
    "  plt.ylabel('Interval Count')\n",
    "  plt.xlabel('time')\n",
    "  plt.show()\n",
    "  plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D2_LinearSystems/solutions/W2D2_Tutorial2_Solution_ecb9dc3a.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 2: Distributional Perspective\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 519
    },
    "colab_type": "code",
    "outputId": "5f098500-bef5-4546-9d9b-03cd5105255b"
   },
   "outputs": [],
   "source": [
    "#@title Video 2: State Transitions\n",
    "# Insert the ID of the corresponding youtube video\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"U6YRhLuRhHg\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "We can run this simulation many times and gather empirical distributions of open/closed states. Alternatively, we can formulate the exact same system probabilistically, keeping track of the probability of being in each state.\n",
    "\n",
    "<!-- Although the system started initially in the Closed ($x=0$) state, over time, it settles into a equilibrium distribution where we can predict on what fraction of time it is Open as a function of the $\\mu$ parameters.  -->\n",
    "\n",
    "(see diagram in lecture)\n",
    "\n",
    "The same system of transitions can then be formulated using a vector of 2 elements as the state vector and a dynamics matrix $\\mathbf{A}$. The result of this formulation is a *state transition matrix*:\n",
    "\n",
    "$\\left[ \\begin{array}{c} C \\\\ O \\end{array} \\right]_{k+1} = \\mathbf{A} \\left[ \\begin{array}{c} C \\\\ O \\end{array} \\right]_k = \\left[ \\begin{array} & 1-\\mu_{\\text{c2o}} & \\mu_{\\text{o2c}} \\\\ \\mu_{\\text{c2o}} & 1-\\mu_{\\text{o2c}} \\end{array} \\right] \\left[ \\begin{array}{c} C \\\\ O \\end{array} \\right]_k$.\n",
    "\n",
    "\n",
    "Each transition probability shown in the matrix is as follows:\n",
    "1. $1-\\mu_{\\text{c2o}}$, the probability that the closed state remains closed. \n",
    "2. $\\mu_{\\text{c2o}}$, the probability that the closed state transitions to the open state.\n",
    "3.  $\\mu_{\\text{o2c}}$, the probability that the open state transitions to the closed state. \n",
    "4. $1-\\mu_{\\text{o2c}}$, the probability that the open state remains open. \n",
    "\n",
    "\n",
    "_Notice_ that this system is written as a discrete step in time, and $\\mathbf{A}$ describes the transition, mapping the state from step $k$ to step $k+1$. This is different from what we did in the exercises above where $\\mathbf{A}$ had described the function from the state to the time derivative of the state.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 2 (2B): Probability Propagation\n",
    "\n",
    "Complete the code below to simulate the propagation of probabilities of closed/open of the ion channel through time. A variable called `x_kp1` (short for, $x$ at timestep $k$ plus 1) should be calculated per each step *k* in the loop. However, you should plot $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "def simulate_prob_prop(A, x0, dt, T):\n",
    "  \"\"\" Simulate the propagation of probabilities given the transition matrix A,\n",
    "  with initial state x0, for a duration of T at timestep dt.\n",
    "\n",
    "  Args:\n",
    "    A (ndarray): state transition matrix\n",
    "    x0 (ndarray): state probabilities at time 0\n",
    "    dt (scalar): timestep of the simulation\n",
    "    T (scalar): total duration of the simulation\n",
    "\n",
    "  Returns:\n",
    "    ndarray, ndarray: `x` for all simulation steps and the time `t` at each step\n",
    "  \"\"\"\n",
    "\n",
    "  # Initialize variables\n",
    "  t = np.arange(0, T, dt)\n",
    "  x = x0 # x at time t_0\n",
    "\n",
    "  # Step through the system in time\n",
    "  for k in range(len(t)-1):\n",
    "      ###################################################################\n",
    "      ## TODO: Insert your code here to compute x_kp1 (x at k plus 1)\n",
    "      raise NotImplementedError(\"Student exercise: need to implement simulation\")\n",
    "      ## hint: use np.dot(a, b) function to compute the dot product\n",
    "      ##       of the transition matrix A and the last state in x\n",
    "      ## hint 2: use np.vstack to append the latest state to x\n",
    "      ###################################################################\n",
    "\n",
    "      # Compute the state of x at time k+1\n",
    "      x_kp1 = ...\n",
    "      # Stack (append) this new state onto x to keep track of x through time steps\n",
    "      x = ...\n",
    "\n",
    "  return x, t\n",
    "\n",
    "# parameters\n",
    "T = 500     # total Time duration\n",
    "dt = 0.1   # timestep of our simulation\n",
    "\n",
    "# same parameters as above\n",
    "# c2o: closed to open rate\n",
    "# o2c: open to closed rate\n",
    "c2o = 0.02\n",
    "o2c = 0.1\n",
    "A = np.array([[1 - c2o*dt, o2c*dt],\n",
    "              [c2o*dt,     1 - o2c*dt]])\n",
    "\n",
    "# initial condition: start as Closed\n",
    "x0 = np.array([[1, 0]])\n",
    "\n",
    "# Uncomment this to plot the probabilities\n",
    "# x, t = simulate_prob_prop(A, x0, dt, T)\n",
    "# plot_state_probabilities(t,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "text",
    "outputId": "b7d16b7a-4130-4b81-d10b-87a035087482"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D2_LinearSystems/solutions/W2D2_Tutorial2_Solution_f426cf5b.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=558 height=414 src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W2D2_LinearSystems/static/W2D2_Tutorial2_Solution_f426cf5b_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Here, we simulated the propagation of probabilities of the ion channel's state changing through time. Using this method is useful in that we can **run the simulation once** and see **how the probabilities propagate throughout time**, rather than re-running and empirically observing the telegraph simulation over and over again. \n",
    "\n",
    "Although the system started initially in the Closed ($x=0$) state, over time, it settles into a equilibrium distribution where we can predict what fraction of time it is Open as a function of the $\\mu$ parameters. We can say that the plot above show this _relaxation towards equilibrium_.\n",
    "\n",
    "Re-calculating our value of the probability of $c2o$ again with this method, we see that this matches the simulation output from the telegraph process! \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 52
    },
    "colab_type": "code",
    "outputId": "08290d22-8447-44d8-ccb3-4c866f1d1fdc"
   },
   "outputs": [],
   "source": [
    "print(\"Probability of state c2o: %.3f\"%(c2o / (c2o + o2c)))\n",
    "x[-1,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 3: Equilibrium of the telegraph process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 519
    },
    "colab_type": "code",
    "outputId": "ff4541b5-56b1-427a-be95-fcd52af8392e"
   },
   "outputs": [],
   "source": [
    "#@title Video 3: Continous vs. Discrete Time Formulation\n",
    "# Insert the ID of the corresponding youtube video\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"csetTTauIh8\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Since we have now modeled the propagation of probabilities by the transition matrix $\\mathbf{A}$ in Section 2, let's connect the behavior of the system at equilibrium with the eigendecomposition of $\\mathbf{A}$.\n",
    "\n",
    "As introduced in the lecture video, the eigenvalues of $\\mathbf{A}$ tell us about the stability of the system, specifically in the directions of the corresponding eigenvectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 69
    },
    "colab_type": "code",
    "outputId": "d29d7e0b-bd1b-425a-9ef2-bb8325006165"
   },
   "outputs": [],
   "source": [
    "# compute the eigendecomposition of A\n",
    "lam, v = np.linalg.eig(A)\n",
    "\n",
    "# print the 2 eigenvalues\n",
    "print(\"Eigenvalues:\",lam)\n",
    "\n",
    "# print the 2 eigenvectors\n",
    "eigenvector1 = v[:,0]\n",
    "eigenvector2 = v[:,1]\n",
    "print(\"Eigenvector 1:\", eigenvector1)\n",
    "print(\"Eigenvector 2:\", eigenvector2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 3 (2C): Finding a stable state\n",
    "\n",
    "Which of these eigenvalues corresponds to the **stable** (equilibrium) solution? What is the eigenvector of this eigenvalue? How does that explain \n",
    "the equilibrium solutions in simulation in Section 2 of this tutorial?\n",
    "\n",
    "_hint_: our simulation is written in terms of probabilities, so they must sum to 1. Therefore, you may also want to rescale the elements of the eigenvector such that they also sum to 1. These can then be directly compared with the probabilities of the states in the simulation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "###################################################################\n",
    "## Insert your thoughts here\n",
    "###################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 52
    },
    "colab_type": "text",
    "outputId": "4dc32baf-8769-4d5b-ec79-aad4e61d086c"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D2_LinearSystems/solutions/W2D2_Tutorial2_Solution_5d5bcbfc.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "In this tutorial, we learned:\n",
    "\n",
    "* The definition of a Markov process with history dependence.\n",
    "* The behavior of a simple 2-state Markov proces--the telegraph process--can be simulated either as a state-change simulation or as a propagation of probability distributions.\n",
    "* The relationship between the stability analysis of a dynamical system expressed either in continuous or discrete time.\n",
    "* The equilibrium behavior of a telegraph process is predictable and can be understood using the same strategy as for deterministic systems in Tutorial 1: by taking the eigendecomposition of the A matrix."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W2D2_Tutorial2",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
